A comparison of the large‐scale gene expression patterns in summer and fall migratory Pantala flavescens (Fabricius) in northern China

Abstract Pantala flavescens (Fabricius) is the most well‐known seasonal migratory insect. This research focused on the molecular response of P. flavescens migration in summer and fall. A total of 17,810 assembled unigenes were obtained and 624 differentially expressed genes (DEGs) were identified in summer migration compared to fall migration. A number of DEGs, including cpr49Ae, itm2b, chitinase, cpr11B, laccase2, nd5, vtg2 and so on, had previously been reported to be involved in cold‐ and high‐temperature resistance. Functional enrichment analysis showed three pathways ‘that antibacterial humoral response, response to bacterial, and lipid transporter activity’ were significantly enriched in summer migration while that six pathways ‘structural constituent of cuticle, chitin binding, mitochondrion, propanoate metabolism, citrate cycle, hypertrophic cardiomyopathy’ were significantly enriched in fall migration. These results will provide a valuable baseline for further understanding of the molecular mechanisms of insect adaptation to different climate migrations.


| INTRODUC TI ON
Insects are ectotherms sensitive to temperature changes.Almost every aspect of insect physiology and behaviour can be influenced by temperature (Cui et al., 2018).Heat and cold stress change insect nervous, respiration, metabolism and endocrine systems.When insects face abiotic environmental stress, such as extreme temperature, they adopt various adaptive strategies, including diapause and migration, to escape unfavourable seasonal conditions (Chapman et al., 2015;Satterfield et al., 2020).
The annual and seasonal migration of insects redistributes nutrients, biomass and species, which has a key influence on essential ecosystem functions (Wotton et al., 2019).Meanwhile, mass-migrating species can interact with other species in places they move through or occupy, which may be important to ecological processes.Every year, billions of dragonflies, butterflies and other insects undertake regular seasonal migrations within and between continents (Wikelski et al., 2006).With a nearly global distribution, Pantala flavescens (Fabricius) (Odonata: Libellulidae) (Figure 1) has the longest documented migratory routine in any known insect and undergoes multigenerational transoceanic migration from India to Africa and back again (Hobson et al., 2012).
Numerous observations by radar studies have revealed that P. flavescens migrate seasonally in northern China, migrating north in summer and south in autumn every year (Cao et al., 2018).Stable isotope tracking of P. flavescens has been conducted to elucidate the origins of migrants (Hobson et al., 2021).In recent years, comparative transcriptomic and next-generation sequencing technologies have been used to explore differences in gene expression between migratory and nonmigratory populations of the species (Talla et al., 2020;Wang et al., 2020;Zhan et al., 2014).For example, some genes related to flight muscle structure, the regulation of hormones and the mobilisation of lipids were revealed in migratory Helicoverpa armigera (Jones et al., 2015).Based on differential expression in flight muscles of migratory and non-migratory monarchs that are correlated with flight metabolic rates, collagen type IV α-1 has been proposed to regulate flight efficiency during long-distance migration (Zhan et al., 2014).The two-way migration of Eastern North American monarchs suggests that epigenetic mechanisms triggered by environmental changes regulate migratory behaviour in this species (Merlin & Liedvogel, 2019).Performing RNA-seq studies in the brains of fall migrants, fall migrants reprogrammed into spring remigrants by cold treatment and wild-caught spring remigrants, could be used to detect differentially expressed genes (DEGs) and non-coding RNAs, as well as cold-dependent RNA splicing or editing events (Merlin et al., 2020).
In contrast to non-migratory P. flavescens, some genes associated with electron transport chain complexes, microtubules, antioxidants and detoxification were found to be highly expressed in migratory P. flavescens (Zhang & Cao, 2023).However, the differences in gene expression profiles in insects between summer and fall migration have rarely been reported.
In this work, our objective was to find a list of the potential DEGs and delineate the function of DEGs between the summer and fall migratory P. flavescens by comparing the gene expression profiles.To further ascertain these different genes, we measured the expression of 10 genes in P. flavescens by quantitative real-time PCR analysis.

| Ethics statement
Pantala flavescens is a migratory insect that is widely distributed across the globe.The insect has been evaluated as 'Least Concern' species by the International Union for Conservation of Nature (IUCN).In China, P. flavescens sampling does not involve endangered or protected species, therefore, no special permits are required for field studies.

| Sample collection, RNA extraction and Illumina sequencing
In July and October 2020, migratory P. flavescens individuals were trapped using searchlight traps on Beihuang Island (38°23′20″ N, 120°54′52″ E), located in the centre of the Bohai Strait in Changdao County, Shandong Province, China (Fu et al., 2014).There is no fresh water at this site, and no larvae of this species were found during daily observations from spring to fall on Beihuang Island.Fifty samples were rinsed three times in 70% ethanol and then stored at −80°C until transported to the laboratory.The samples were divided into summer migratory individuals trapped in July (M7) and fall migratory individuals trapped in October (M10).
Total RNA was extracted from the thoracic muscle of each sample using TRIzol® Reagent (Invitrogen) according to the manufacturer's instructions.The purity, concentration and integrity of the RNA were verified by a Nanodrop2000 Spectrophotometer (Thermo Fisher Scientific) and 1% agarose gel electrophoresis.Samples with three biological replicates from summer and fall migratory individuals were labelled M7-1, M7-2 and M7-3 and M10-1, M10-2 and M10-3.The RNA used in each biological repetition comes from a mixed RNA of five individuals.In total, six mixed RNA samples were sent to Shanghai Majorbio Bio-Pharm Technology Co., Ltd.(Shanghai, China) for cDNA library construction and RNA-seq.

| De novo assembly and gene annotation
The raw Illumina sequence data were trimmed and controlled for quality by fastp (https:// github.com/ OpenG ene/ fastp ) with default parameters (Chen et al., 2018).After filtering, the remaining reads were deemed clean reads and stored in FASTQ format.The clean reads were assembled by using the Trinity assembly program with default parameters (Grabherr et al., 2011).TGICL clustering software v2.1 with default parameters was used to cluster and remove redundancy to obtain nonredundant unigenes (Pertea et al., 2003).The quality of the unigenes was measured using TransRate (v2.4.0).The integrity of the transcriptome assembly was assessed by BUSCO v3.0.2 (Simão et al., 2015).
To annotate the dragonfly transcriptome, we performed a BLAST v2.2.23 search with a cut-off E-value of 10 −5 against six databases: the nonredundant protein sequence (nr) and nucleotide sequence (nt) databases, the KEGG database, the Gene Ontology (GO) database, a manually annotated and reviewed protein sequence database (Swiss-Prot) and the COG database.The sequence orientations and coding sequences of the unigenes were determined based on the BLASTx best hits from these databases.Gene Ontology annotation of unigenes was conducted by Blast2GO v2.5.0 software through a search of the nr database (Conesa et al., 2005), and WEGO software was used to determine GO functional classification and assess the functional distribution (Ye et al., 2006).

F I G U R E 1
The photo of P. flavescens.

| Abundance estimation and identification of DEGs
The unigenes from each sample were compared with the reference unigene database using Bowtie2 v2.1.0software (Langmead & Salzberg, 2012).The expression level of the unigenes was estimated using RSEM v1.2.12 (Li & Dewey, 2011) with the default parameters for each sample, employing the reads per kilo base per million mapped reads (RPKM) values.DESeq2 was applied for the identification of DEGs in different samples (Love et al., 2014).Through comparison, we used the absolute value of log 2 ratio ≥ 2 and p-adjust ≤.01 as the threshold to judge the significance level of gene expression differences.GO functional enrichment and KEGG pathway enrichment analyses were carried out by Goatools (https:// github.com/ tangh aibao/ Goatools) and KOBAS (http:// kobas.cbi.pku.edu.cn) (Xie et al., 2011).Gene clustering analysis was performed to cluster DEGs with RSEM with p-adjust ≤.05 (Li & Dewey, 2011).

| Quantitative real-time PCR analysis (qRT-PCR validation)
To validate the reliability of the RNA-Seq data, we randomly selected 10 enriched DEGs for qRT-PCR analysis with β-actin as the reference gene (Ridgeway & Timm, 2015;Zhao et al., 2022).Primers for qRT-PCR are shown in Table 1.The total RNA of each sample was extracted using TRIzol® Reagent (Invitrogen) according to the manufacturer's protocol.HiScript Q RT SuperMix for qPCR (Vazyme, China) was used to synthesise cDNA.qRT-PCR was carried out in a LightCycler® 480 Real-time PCR system (Roche) using ChamQ Universal SYBR qPCR Master Mix (Vazyme, China).Each reaction was conducted with three independent biological replicates.The relative expression ratio of the target gene was determined using Ct (threshold cycles) and calculated using the 2 −∆∆Ct method (Livak & Schmittgen, 2001).The qRT-PCR data were analysed using SPSS v16.0 software.

| Sequencing and de novo assembly
After filtering low-quality sequences, the clean data for each sample were more than 7.0 Gb, with a total of 316.79 million clean reads, and 46.62 Gb were obtained with 98.3% Q20 and 94.74% Q30 (Table 2).Through the assembly and alignment of the clean reads, 17,810 unigenes and 27,701 transcripts were obtained with an average length of 2046 bp and an N50 value of 3583 bp.
For the 17,810 assembled unigenes, a total of 10,920 unigenes (61.32%) were assigned annotations in at least one database, while the remaining unigenes (38.68%) were of undetermined function.in the GO database, 7189 unigenes (40.37%) to genes in the KEGG database, 10,076 unigenes (56.58%) to genes in the eggNOG database, 10,510 unigenes (59.02%) to genes in the NR database, 9050 unigenes (50.82%) to genes in the Swiss-Prot database and 9496 unigenes (53.32%) to genes in the Pfam database (Table 3).

| Gene expression differences between summer and fall migrants
A total of 624 DEGs were identified by DESeq2 in P. flavescens between summer and fall migration with the criterion of log 2 ratio ≥ 2 and a false discovery rate (FDR) of ≤0.01.Among 624 DEGs, there were 352 significantly upregulated and 272 downregulated DEGs in the M7 group compared to the M10 group (Figure 2).The top 15 up-and downregulated genes were chosen by the highest expression changes FC (M10/ M7) (Table 4).The putative functions of these genes were identified according to the NR database.Among the 15 upregulated genes, the proved that the RNA-seq data were reliable (Figure 3).

| Functional gene enrichment analysis
In vivo, many genes often form networks that coordinate with each other to perform biological functions.Enrichment pathway analysis is helpful to understand the biological function of genes.GO enrichment analysis revealed that 10 GO terms were significantly enriched in the M7 group versus the M10 group, among which the molecular functions (MF) were associated with structural constituent of cuticle, chitin binding, lipid transporter activity and iron-ion binding; the cell components (CC) were the mitochondria and extracellular region; and the biological processes (BP) were response to bacterium and humoral response (Table 5).
KEGG is a database resource for the utility of biological systems, such as cells, organisms and ecosystems, and has been constructed from genomic and molecular-level information.In this study, KEGG enrichment analyses showed that four pathways were significantly enriched, including hypertrophic cardiomyopathy, propanoate metabolism, citrate cycle and dilated cardiomyopathy (Table 5).

| Functional gene cluster analysis
To understand the expression level of different genes in the enrichment pathway, we performed cluster analysis of these genes.The  genes involved in the structural constituent of cuticle and chitin were overexpressed, while the expression level of Chitinase was downregulated in fall migrants (Figure 4).In the four categories 'response to bacterium, defence response to bacterium, antibacterial/antimicrobial humoral response', the genes annotated for sarcotoxin-2A and defensins were overexpressed in summer migrants (Figure 5a).For the GO term 'lipid transporter activity', the expression levels of vtg2encoded vitellogenin and fasn2-encoded fatty acid synthase 2 were upregulated in summer migrants (Figure 5b).For the GO term 'mitochondrion', the genes associated with mitochondrial protein synthesis and energy production were overexpressed in fall migrants (Figure 6a).
In the propanoate metabolism and citrate cycle enrichment pathways, the expression levels of the genes related to energy production were upregulated significantly in fall migrants (Figure 7).

| DISCUSS ION
The decision to migrate is simultaneously affected by external environmental conditions and the insects' own behaviour and physiological state (Chapman et al., 2015).Seasonal ecosystem linkage was a key ecosystem property for maintaining migratory polymorphism in partially migratory animals (Tanaka et al., 2023).In this study, we identified 624 DEGs, 10 enriched GO pathways and 4 enriched KEEG pathways in the M7 group versus the M10 group.These enrichment pathways are mainly involved in cellular exoskeleton composition, bacterium and humoral response, energy metabolism and muscle contraction.Our findings will help reveal new genes associated with the transoceanic migration of other insects and their impact on global ecology.We have ordered our discussion below based on the DEGs in enriched pathways.
TA B L E 5 List of the significant results from the gene ontology enrichment analysis and KEGG enrichment analysis between M7 and M10.

| DEGs related to exoskeleton composition
In our study, the expression level of chitinase was significantly upregulated in summer migrants (Figure 4a), while the genes annotated for cuticle protein, structural constituent of cuticle and endocuticle structural glycoproteins were significantly upregulated in fall migrants (Figure 4b).Chitin is one of the major components of insect cuticle, which is already known to be a dynamic structure, involved in protecting against water loss, chemical erosion, physical damage and pathogen invasion (Muthukrishnan et al., 2020;Qu et al., 2017).The cuticular structure-associated transcripts were enriched in many insects exposed to low temperature, suggesting that insect cuticle may play an important role in adaptation to low temperature (Dunning et al., 2014;Huang et al., 2017;Zhou et al., 2019).The individuals in the M10 group

| DEGs related to defence and detoxification
Defensins and sarcotoxins are both types of antibacterial proteins and have a large spectrum of antimicrobial action.Defensins are crucial components of the immune response of insects against microbial invasion, parasitic and pathogenic infections (Pasupuleti et al., 2012;Raj & Dentino, 2002).The expression of sarcotoxin II increased strongly in fat bodies after injury and it has been shown to be effective only against a few gram-negative bacteria (Andoh et al., 2018).
The expression of AMP genes including def1-encoded defensin and loc108905441-encoded sarcotoxin II, which protect against bacteria and harmful microorganisms and repair environmental damage, was overexpressed in M7 (Figure 5a), which may be related to the loose binding of the exoskeleton of individuals in this group.
Cytochrome P450 (P450s) are the main detoxifying enzymes involved mainly in the metabolism of a variety of endogenous and exogenous compounds and agrochemical insecticide resistance in insects (Feyereisen, 2012;Jeschke, 2016).For example, overexpression of CYP 6G1 in insects resulted in the increased metabolic detoxification of insecticides (Huang et al., 2013;Joussen et al., 2008).
The upregulation of cyp expression in M7 group may be related to their higher exposure to toxic substances (Figure 5a).

| DEGs related to energy consumption
Vitellogenin encodes the major yolk protein precursor and plays an essential role in the reproduction and spread of all oviparous species, including insects (Tufail et al., 2010).Vitellogenin is not only important for bee reproduction but also involved in their immune response (Harwood et al., 2019;Sun & Zhang, 2015), antioxidative stress response (Seehuus et al., 2006;Salmela et al., 2016) and metal ion transport (Amdam et al., 2004).In the lipid transport term, the expression level of vtg2 was greatly upregulated in the M7 group, which may have been related to storage of energy, antioxidant activity and metal ion transport in the individuals of this group (Figure 5b).The expression levels of fasn2-and apoe-encoded apolipoprotein were upregulated in the M7 group, which may indicate that individuals stored more fat in this group (Figure 5b).
Mitochondria are essential organelles that execute and coordinate various metabolic pathways involved in ATP synthesis through the tricarboxylic acid cycle and oxidative phosphorylation in the cell (Ng et al., 2021).The genes involved in ribosome assembly and protein synthesis were overexpressed in M10, which indicates the synthesis of a great number of proteins and ATP production in the mitochondria (Figure 6a).Trifunctional enzyme subunits α and β compose mitochondrial trifunctional proteins that catalyse the last three steps in long-chain fatty acid metabolism for energy supply (van Vliet et al., 2018).The expression levels of mtpbeta and acads encoded short-chain acyl-CoA dehydrogenase upregulated in the M10 group can promote fatty acid metabolism (Figure 6b).The pyruvate dehydrogenase complex composed of pyruvate dehydrogenase, dihydrolipoyl transacetylase and dihydrolipoyl dehydrogenase catalyses, the oxidative decarboxylation of pyruvate to acetyl coenzyme A, linking glycolysis to the tricarboxylic acid cycle (Holness & Sugden, 2003).The expression levels of pdha and dld-1 were significantly upregulated in the M10 group, which can accelerate glycolysis, the tricarboxylic acid cycle and ATP production (Figure 6c).

| DEGs related to muscle contraction
Sarcomeric myosin heavy-chain (MYH) genes are expressed in the major cardiac and skeletal muscles.Myosin produces the force necessary for a variety of cellular movements by hydrolysing ATP and interacting with actin filaments (Lee et al., 2019).The high expression levels of mhc and actin indicated that contractions and cellular movements of individuals in the M10 group were more frequent and intense (Figure 7).The cardiac isoform of the sarcoplasmic/ F I G U R E 7 Comparative expression levels of DEGs associated with dilated cardiomyopathy.The colour scale, which encompasses the lowest (blue) to the highest (red) TPM value, is displayed on the right side.
vast majority of upregulated genes were associated with structural constituent of cuticle (11/15).The 15 top downregulated genes were mainly involved in vitellogenin and lipid transporter activity (6/15), defensin and sarcotoxin (3/15) and ion binding (3/15) (Table 4).To validate the RNA-Seq data, quantitative PCR (qPCR) was used to measure the expression levels of 10 DEGs.The results of the qRT-PCR validation

F I G U R E 2
Volcano plots of DEGs in the M7 group versus M10.The x axis represents the log 2 -fold change, and the y axis represents the negative log 10 of the corrected p value.Blue points denote significantly downregulated genes, and red denotes significantly upregulated genes in each comparison with the criteria of |log 2 (Foldchange)| > 2 and corrected p value <.01.
and conversion, such as mtpbeta-encoded trifunctional enzyme subunit alpha, dld-1-encoded dihydrolipoyl dehydrogenase, hibch-encoded 3-hydroxyisobutyry1-COA hydrolase and pdha-encoded pyruvate dehydrogenase, were upregulated in fall migrants (Figure 6b,c).These results indicated that the function of mitochondria was enhanced to release energy when P. flavescens were under cold stress.In the hypertrophic cardiomyopathy and dilated cardiomyopathy enrichment pathways, the genes related to movement, such as mhc-encoded myosin heavy chain, actin, and mca-3-encoded calcium-transporting ATPase,

F
I G U R E 3 qRT-PCR validation of DEGs.Three replicates were analysed.The X-axis represents different unigenes.The selected genes are listed in Table1.The Y-axis shows the relative expression levels of genes.Actin was used as a reference gene.The error bars indicate standard deviations.

F
Comparative expression levels of DEGs related to chitin (a) and cuticle (b) in M7 versus M10.The colour scale, which encompasses the lowest (blue) to the highest (red) TPM value, is displayed on the right side.

F
Comparative expression levels of DEGs related to defence (a) and lipid transport (b) in M7 versus M10.The colour scale, which encompasses the lowest (blue) to the highest (red) TPM value, is displayed on the right side.F I G U R E 6 Comparative expression levels of DEGs associated with mitochondrial metabolism (a), propanoate metabolism (b) and citrate cycle (c) in the M7 group versus M10.The colour scale, which encompasses the lowest (blue) to the highest (red) TPM value, is displayed on the right side.had more compact chitin and cuticle layers, protecting the dragonflies against freezing and harm from the external environment.The upregulated expression of chitinase can destroy the compact structure of outer epidermis and cuticle of individuals in the M7 group, and make its surface loose to heat dissipation and water evaporation which is suitable for the high temperature in July.
Summary of the transcriptome data filtering.Summary statistics of annotation of all unigenes.
TA B L E 2 Top 15 upregulated and downregulated genes.List are ordered by FC (M10/M7).
TA B L E 4